Libraries:
# libraries -------------------
library(dplyr)
library(tidyverse)
library(ggplot2)
library(GGally)
library(glmmTMB)
library(TMB)
library(emmeans)
library(DHARMa)
library(lmtest)
Set working directory
knitr::opts_knit$set(root.dir = '/Users/user/Desktop/Data/Regen_RProj/')
Functions:
# functions -------------------
Source large seedling data:
source("Scripts/LS_Import.R")
Pivot wider to create dataframe where each row is for one plot and has total details for each species group
ls_merge2 <- ls_merge %>%
select(Plot_No, Region, Treat_Type, Site, Species_Groups, Total) #this is dropping browse and stump sprout data
ls_merge2 <- ls_merge2 %>%
pivot_wider(names_from = Species_Groups, values_from = Total)
Import time since data and add it to the large seedling dataset
time_since <- read_csv("CleanData/Treat_Year_Data.csv")
ls_merge3 <- merge(ls_merge2, time_since, by = 'Site')
#log transform time from treatment data
ls_merge3$l.TFT <- log(ls_merge3$Time_from_Treat)
Run the ‘Add_BA’ script and merge with dataset:
source("Scripts/Add_BA.R")
# merge with ls dataset -------------------
ls_merge4 <- merge(ls_merge3, prism_BA, by = "Plot_No")
Run ‘Ground_Data.R’ script and add it to large seedling dataset:
source("Scripts/Ground_Data.R")
# merge with ls dataset -------------------
ls_merge5 <- merge(ls_merge4, ground3, by = "Plot_No")
Source and add veg data
source("Scripts/Veg_Data.R")
# merge with ls dataset
ls.all <- merge(ls_merge5, veg3, by = "Plot_No")
rm(ls_merge5,
ls_merge2,
ls_merge3,
ls_merge4)
The large seedling count data is taken in 10m2 plots; basal area is measured in hectares; veg and soil data is taken in 1m2 plots.
Large seedling data will be converted into 1m2 plots in order to compare across and reduce the amount of scales of data collection to two: 1m2 plots and per hectare observations.
ls.all$PIRI.1m <- ls.all$PIRI/10
ls.all$SO.1m <- ls.all$Shrub_Oak/10
ls.all$Other.1m <- ls.all$Other/10
Create log transformed categories for newly added variables, then select for just the desired variables:
ls.all$l.PIRI1 <- log(ls.all$PIRI.1m + 1)
ls.all$l.SO1 <- log(ls.all$SO.1m + 1)
ls.all$l.other1 <- log(ls.all$Other.1m + 1)
ls.all2 <- ls.all %>%
select(Treat_Type, Region, Site, Plot_No, PIRI, PIRI.1m, l.PIRI1, Shrub_Oak, SO.1m, l.SO1, Other, Other.1m, l.other1, Time_from_Treat, l.TFT, BA_HA, l.BA_HA, PIRI.BA_HA, l.BA_piri, Mineral_Soil, l.Mineral, Litter_Duff, avgLD, avgLD_l, Veg_Total, l.Veg_Total) %>%
arrange(Treat_Type)
Select just for numerical vs log and then look at paired plots:
#not transformed
ls.num <- ls.all2 %>%
select(PIRI, Shrub_Oak, Other, Time_from_Treat, BA_HA, PIRI.BA_HA, Mineral_Soil, avgLD, Veg_Total, Treat_Type)
ggpairs(ls.num)
ggpairs(ls.num, aes(color = Treat_Type))
#log transformed
ls.numl <- ls.all2 %>%
select(l.PIRI, l.SO, l.other, l.TFT, l.BA_HA, l.BA_piri, l.Mineral, avgLD_l, l.Veg_Total, Treat_Type)
ggpairs(ls.numl)
ggpairs(ls.numl, aes(color = Treat_Type))
rm(ls.num,
ls.numl)
Can see the correlation coefficients for linear (Pearsons) relationships. None of them appear very strong, except for ones that are analogs (avg LD vs mineral soil exposure; ba/ha vs piri ba/ha)
Log transformed average litter depth and basal area per hectare have a weak relationship (corr 0.33), which does make sense.
Full dataset is called ls.all2
Starting first with model analysis without treatment type:
ls.m1 <- glmmTMB(PIRI.1m ~ l.SO1 + l.other1 + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all2,
family = poisson)
## Warning in glmmTMB(PIRI.1m ~ l.SO1 + l.other1 + l.BA_HA + l.BA_piri + l.Mineral
## + : non-integer counts in a poisson model
AIC(ls.m1) #98.5
## [1] 98.45309
# gives warning about non-integer values (as i've reduced the counts (/10) to get them to the 1m2 ) - if BA isn't important, I could try again with no BA in dataset and keep 10m2 observations ....
summary(ls.m1)
## Family: poisson ( log )
## Formula:
## PIRI.1m ~ l.SO1 + l.other1 + l.BA_HA + l.BA_piri + l.Mineral +
## avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Data: ls.all2
##
## AIC BIC logLik deviance df.resid
## 98.5 140.6 -39.2 78.5 488
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot_No:Site (Intercept) 2.561e-09 5.061e-05
## Site (Intercept) 6.410e+00 2.532e+00
## Number of obs: 498, groups: Plot_No:Site, 498; Site, 47
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.69499 3.33513 -1.408 0.159
## l.SO1 -1.13705 1.00373 -1.133 0.257
## l.other1 -0.21927 1.53411 -0.143 0.886
## l.BA_HA -0.28410 0.59909 -0.474 0.635
## l.BA_piri -0.38588 0.60304 -0.640 0.522
## l.Mineral 0.08863 0.34845 0.254 0.799
## avgLD_l -0.08137 0.87393 -0.093 0.926
## l.Veg_Total -0.33215 0.65895 -0.504 0.614
ls.m2 <- glmmTMB(PIRI.1m ~ l.SO1 + l.other1 + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all2,
family = poisson)
## Warning in glmmTMB(PIRI.1m ~ l.SO1 + l.other1 + l.BA_HA + l.Mineral + avgLD_l +
## : non-integer counts in a poisson model
AIC(ls.m2) #96.8
## [1] 96.83059
lrtest(ls.m1, ls.m2) #p = 0.5
## Likelihood ratio test
##
## Model 1: PIRI.1m ~ l.SO1 + l.other1 + l.BA_HA + l.BA_piri + l.Mineral +
## avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI.1m ~ l.SO1 + l.other1 + l.BA_HA + l.Mineral + avgLD_l +
## l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 10 -39.227
## 2 9 -39.415 -1 0.3775 0.5389
ls.m3 <- glmmTMB(PIRI.1m ~ l.SO1 + l.other1 + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all2,
family = poisson)
## Warning in glmmTMB(PIRI.1m ~ l.SO1 + l.other1 + l.Mineral + avgLD_l +
## l.Veg_Total + : non-integer counts in a poisson model
AIC(ls.m3) #96.0
## [1] 96.00496
lrtest(ls.m2, ls.m3) # p = 0.3
## Likelihood ratio test
##
## Model 1: PIRI.1m ~ l.SO1 + l.other1 + l.BA_HA + l.Mineral + avgLD_l +
## l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI.1m ~ l.SO1 + l.other1 + l.Mineral + avgLD_l + l.Veg_Total +
## offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 9 -39.415
## 2 8 -40.002 -1 1.1744 0.2785
rm(ls.m1, ls.m2, ls.m3)
ls.m11a <- glmmTMB(PIRI.1m ~ avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all2,
family = poisson)
## Warning in glmmTMB(PIRI.1m ~ avgLD_l + offset(l.TFT) + (1 | Site/Plot_No), :
## non-integer counts in a poisson model
AIC(ls.m11a) #AIC = 89.3, was just interested to see the AIC of this model vs. ls at 10m2 observations
## [1] 89.37025
#i'm just interested to see model fit on this 1m model
ls.m11a_sr <- simulateResiduals(ls.m11a, n = 1000, plot = T) #fails very much so
Basal area doesn’t seem important, so I’m going to go back to the 10m2 and 1m2 scales and lose per hectare observations. I’ll need to rework the dataset
Revised data set with LS observations at 10m2 scale; this means no non-interger values for the poisson distro
ls.all3 <- ls.all
ls.all3$l.PIRI <- log(ls.all3$PIRI + 1)
ls.all3$l.SO <- log(ls.all3$Shrub_Oak + 1)
ls.all3$l.other <- log(ls.all3$Other + 1)
ls.all3 <- ls.all3 %>%
select(Treat_Type, Region, Site, Plot_No, PIRI, l.PIRI, Shrub_Oak, l.SO, Other, l.other, Time_from_Treat, l.TFT, BA_HA, l.BA_HA, PIRI.BA_HA, l.BA_piri, Mineral_Soil, l.Mineral, Litter_Duff, avgLD, avgLD_l, Veg_Total, l.Veg_Total) %>%
arrange(Treat_Type)
Begin modeling again:
#double checking about basal area variables, even though this is at 3 scales
ls.m4 <- glmmTMB(PIRI ~ l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(ls.m4) #AIC is 291
#test piri ba
ls.m5 <- glmmTMB(PIRI ~ l.SO + l.other + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(ls.m5) #289.2
lrtest(ls.m4, ls.m5) #p = 0.7
#test ba
ls.m6 <- glmmTMB(PIRI ~ l.SO + l.other + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(ls.m6) #287.3
lrtest(ls.m5, ls.m6) #p=0.8
#test mineral soil
ls.m7 <- glmmTMB(PIRI ~ l.SO + l.other + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(ls.m7) #285.3
lrtest(ls.m6, ls.m7) # p = 0.8
#test litter depth
ls.m8 <- glmmTMB(PIRI ~ l.SO + l.other + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(ls.m8) #289.5
lrtest(ls.m7, ls.m8) # p = 0.01, keep avg LD
#return avg ld, test other
ls.m9 <- glmmTMB(PIRI ~ l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(ls.m9) #285.1
lrtest(ls.m7, ls.m9) #p = 0.2
# test shrub oak
ls.m10 <- glmmTMB(PIRI ~ avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(ls.m10) #285.2
lrtest(ls.m10, ls.m9) #p = 0.1
#test veg cover ------ seems like this is the best model
ls.m11 <- glmmTMB(PIRI ~ avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(ls.m11) # p = 283.2
lrtest(ls.m10, ls.m11) # p = 0.97
rm(ls.m4, ls.m5, ls.m6, ls.m7, ls.m8, ls.m9, ls.m10)
Now to test model fit:
ls.m11_sr <- simulateResiduals(ls.m11, n = 1000, plot = TRUE) #doesn't pass ks
testResiduals(ls.m11_sr)
testZeroInflation(ls.m11_sr)
# I wonder if I added treat type, if the model would pass
ls.m12 <- glmmTMB(PIRI ~ Treat_Type + avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(ls.m12) #284.6
lrtest(ls.m12, ls.m11) #says treat type isn't important
ls.m12_sr <- simulateResiduals(ls.m12, n = 1000, plot = T) #looks better, still not perfect
testResiduals(ls.m12_sr) #passes
testZeroInflation(ls.m12_sr) #passes
testQuantiles(ls.m12_sr)
#going to test a model with more variables (and not treat type), to see if that fits better
ls.m10_sr <- simulateResiduals(ls.m10, n = 1000, plot = T)
#tests on model 9 failed & on 10
Models without treatment type failed in model fit. Going to run variable elimination again, starting with models with treatment type
ls.m13 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(ls.m13) #293.2
## [1] 293.1687
ls.m14 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(ls.m14) #291.3
## [1] 291.3429
ls.m15 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(ls.m15) #289.5
## [1] 289.4573
ls.m16 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(ls.m16) #287.5
## [1] 287.4731
ls.m17 <- glmmTMB(PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(ls.m17) #286.3
## [1] 286.315
lrtest(ls.m13, ls.m14, ls.m15, ls.m16, ls.m17)
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral +
## avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.Mineral + avgLD_l +
## l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 3: PIRI ~ Treat_Type + l.SO + l.other + l.Mineral + avgLD_l + l.Veg_Total +
## offset(l.TFT) + (1 | Site/Plot_No)
## Model 4: PIRI ~ Treat_Type + l.SO + l.other + avgLD_l + l.Veg_Total +
## offset(l.TFT) + (1 | Site/Plot_No)
## Model 5: PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) +
## (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 14 -132.58
## 2 13 -132.67 -1 0.1742 0.6764
## 3 12 -132.73 -1 0.1144 0.7352
## 4 11 -132.74 -1 0.0159 0.8998
## 5 10 -133.16 -1 0.8418 0.3589
ls.m18 <- glmmTMB(PIRI ~ Treat_Type + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(ls.m18) #286.6
## [1] 286.5651
lrtest(ls.m17, ls.m18) #p = 0.13
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) +
## (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + avgLD_l + l.Veg_Total + offset(l.TFT) + (1 |
## Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 10 -133.16
## 2 9 -134.28 -1 2.2502 0.1336
ls.m19 <- glmmTMB(PIRI ~ Treat_Type + avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
summary(ls.m19) #284.6
## Family: poisson ( log )
## Formula:
## PIRI ~ Treat_Type + avgLD_l + offset(l.TFT) + (1 | Site/Plot_No)
## Data: ls.all3
##
## AIC BIC logLik deviance df.resid
## 284.6 318.3 -134.3 268.6 490
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot_No:Site (Intercept) 2.096 1.448
## Site (Intercept) 9.911 3.148
## Number of obs: 498, groups: Plot_No:Site, 498; Site, 47
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.7216 1.8797 -4.640 3.49e-06 ***
## Treat_TypeFallRx 2.8500 1.8888 1.509 0.13133
## Treat_TypeHarvest 5.7603 2.0508 2.809 0.00497 **
## Treat_TypeMowRx 2.1656 1.8477 1.172 0.24117
## Treat_TypeSpringRx 2.7379 2.1166 1.294 0.19582
## avgLD_l -1.1665 0.5173 -2.255 0.02413 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lrtest(ls.m18, ls.m19) #0.98
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + avgLD_l + l.Veg_Total + offset(l.TFT) + (1 |
## Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + avgLD_l + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 9 -134.28
## 2 8 -134.28 -1 6e-04 0.9812
ls.m20 <- glmmTMB(PIRI ~ Treat_Type + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(ls.m20) #288
## [1] 288.0281
lrtest(ls.m19, ls.m20) # p = 0.02
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + avgLD_l + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 8 -134.28
## 2 7 -137.01 -1 5.4624 0.01943 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rm(ls.m13, ls.m14, ls.m15, ls.m16, ls.m17, ls.m18, ls.m20)
Ok, test model fit with treat type (same as model 11)
ls.m19_sr <- simulateResiduals(ls.m19, n = 1000, plot = TRUE) #passes KS, quantile deviations fails, but it still could be an accepted model
testResiduals(ls.m19_sr) #pases
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.055868, p-value = 0.0893
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.6209e-05, p-value = 0.556
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.009086345
## sample estimates:
## outlier frequency (expected: 0.00154618473895582 )
## 0
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.055868, p-value = 0.0893
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.6209e-05, p-value = 0.556
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.009086345
## sample estimates:
## outlier frequency (expected: 0.00154618473895582 )
## 0
testZeroInflation(ls.m19_sr) #passes
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.005, p-value = 0.94
## alternative hypothesis: two.sided
testQuantiles(ls.m19_sr) #fails
##
## Test for location of quantiles via qgam
##
## data: simulationOutput
## p-value = 0.00256
## alternative hypothesis: both
Trying pairwise comparison of model fit:
emmeans(ls.m19, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
## $emmeans
## Treat_Type rate SE df asymp.LCL asymp.UCL
## Control 0.000168 0.000292 Inf 5.63e-06 0.00503
## FallRx 0.002912 0.004627 Inf 1.29e-04 0.06559
## Harvest 0.053471 0.089469 Inf 2.01e-03 1.42029
## MowRx 0.001469 0.002266 Inf 7.14e-05 0.03022
## SpringRx 0.002603 0.004432 Inf 9.25e-05 0.07322
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## Control / FallRx 0.05785 0.10926 Inf 1 -1.509 0.5566
## Control / Harvest 0.00315 0.00646 Inf 1 -2.809 0.0399
## Control / MowRx 0.11468 0.21190 Inf 1 -1.172 0.7674
## Control / SpringRx 0.06471 0.13696 Inf 1 -1.294 0.6952
## FallRx / Harvest 0.05446 0.10554 Inf 1 -1.502 0.5613
## FallRx / MowRx 1.98261 3.42698 Inf 1 0.396 0.9948
## FallRx / SpringRx 1.11861 2.24433 Inf 1 0.056 1.0000
## Harvest / MowRx 36.40752 69.35922 Inf 1 1.887 0.3243
## Harvest / SpringRx 20.54164 43.83726 Inf 1 1.416 0.6171
## MowRx / SpringRx 0.56421 1.10087 Inf 1 -0.293 0.9984
##
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log scale
Control vs Harvest is the only significant different (p = 0.0399)
I’d like to see some graphs before I close out
library(sjPlot)
library(sjlabelled)
##
## Attaching package: 'sjlabelled'
## The following object is masked from 'package:forcats':
##
## as_factor
## The following object is masked from 'package:ggplot2':
##
## as_label
## The following object is masked from 'package:dplyr':
##
## as_label
library(sjmisc)
##
## Attaching package: 'sjmisc'
## The following object is masked from 'package:purrr':
##
## is_empty
## The following object is masked from 'package:tidyr':
##
## replace_na
## The following object is masked from 'package:tibble':
##
## add_case
set_theme(base = theme_classic(),
theme.font = 'serif',
axis.title.size = 1.5,
axis.textsize.x = 1.5,
axis.textsize.y = 1.5,
title.size = 2.5,
title.align = "center",
legend.pos = "right",
legend.size = 1.5,
legend.title.size = 1.5,
#legend.bordercol = "black",
legend.item.size = .75)
# plot_model(ls.m19) #this is incidence rate ratios
#
# plot_model(ls.m19, type = "diag") #random vs normal quantiles
#
# plot_model(ls.m19, type = "re") #this plots random effects
#
# plot_model(ls.m19,
# type = 'pred',
# terms = 'Treat_Type') #plot marginal effects (i might have done this wrong)
# LS PIRI plot 1 ------------------- LD vs. TT
plot_model(ls.m19,
type = 'pred',
terms = c('avgLD_l', 'Treat_Type'),
axis.title = c("Average Leaf Litter Depth (log transformed)", "Total Count of Pitch Pine"),
title = "Predicted Counts of Pitch Pine Seedlings \n >/=50cm & <2.5cm DBH",
legend.title = "Treatment Type",
line.size = 1,
value.offset = 'Treat_Type',
ci.lvl = NA,
colors = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))
## Warning in checkTerms(data.tmb1$terms, data.tmb0$terms): Predicting new random effect levels for terms: 1 | Plot_No:Site
## Disable this warning with 'allow.new.levels=TRUE'
# LS PIRI plot 2 ------------------- Time from Treatment vs TT
plot_model(ls.m19,
type = 'pred',
terms = c('l.TFT', 'Treat_Type'),
axis.title = c("Time from Treatment (log transformed)", "Total Count of Pitch Pine"),
title = "Predicted Counts of Pitch Pine Seedlings \n >/=50cm & <2.5cm DBH",
legend.title = "Treatment Type",
line.size = 1,
value.offset = 'Treat_Type',
ci.lvl = NA,
colors = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))
## Warning in checkTerms(data.tmb1$terms, data.tmb0$terms): Predicting new random effect levels for terms: 1 | Plot_No:Site
## Disable this warning with 'allow.new.levels=TRUE'
Should I think about graphs where the offset is included?
Random slope for models?
so.ls1 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(so.ls1) # 3377.9
## [1] 3377.913
# test piri ba
so.ls2 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(so.ls2) #3379.8
## [1] 3379.761
lrtest(so.ls1, so.ls2) #p = 0.0498, so maybe
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.BA_piri +
## l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1 |
## Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral +
## avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 14 -1675.0
## 2 13 -1676.9 -1 3.8484 0.04979 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# test total BA
so.ls3 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(so.ls3) #3379.8
## [1] 3379.829
lrtest(so.ls1, so.ls3) # p = 0.0478, keep BA for now
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.BA_piri +
## l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1 |
## Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_piri + l.Mineral +
## avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 14 -1675.0
## 2 13 -1676.9 -1 3.916 0.04783 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# test mineral
so.ls4 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(so.ls4) #3381.7
## [1] 3381.666
lrtest(so.ls2, so.ls4) # p = 0.048
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral +
## avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + avgLD_l +
## l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 13 -1676.9
## 2 12 -1678.8 -1 3.905 0.04814 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# test avg ld
so.ls5 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(so.ls5) #3378.3 ---- keep mineral, lower AIC
## [1] 3378.265
lrtest(so.ls5, so.ls2) # p = 0.48
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral +
## l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral +
## avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 12 -1677.1
## 2 13 -1676.9 1 0.504 0.4777
# test piri
so.ls6 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(so.ls6) #3380
## [1] 3380.052
lrtest(so.ls6, so.ls5) # p = 0.051
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total +
## offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral +
## l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 11 -1679.0
## 2 12 -1677.1 1 3.7868 0.05166 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# test other
so.ls7 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(so.ls7) #3388.2
## [1] 3388.225
lrtest(so.ls6, so.ls7) # p = 0.001
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total +
## offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.BA_HA + l.Mineral + l.Veg_Total +
## offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 11 -1679.0
## 2 10 -1684.1 -1 10.173 0.001425 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# test veg
so.ls8 <- glmmTMB(Shrub_Oak ~ Treat_Type +l.other + l.BA_HA + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(so.ls8) # 3383.4
## [1] 3393.376
lrtest(so.ls8, so.ls6) # p > 0.001
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + offset(l.TFT) +
## (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total +
## offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 10 -1686.7
## 2 11 -1679.0 1 15.324 9.055e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Seemingly best models
so.ls5 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(so.ls5) #3378.3 ---- keep mineral, lower AIC
## [1] 3378.265
so.ls5_sr <- simulateResiduals(so.ls5, n = 1000, plot = TRUE) #quantile test not looking good
testResiduals(so.ls5_sr)
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.055614, p-value = 0.09186
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.0054661, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.00813253
## sample estimates:
## outlier frequency (expected: 0.0015863453815261 )
## 0
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.055614, p-value = 0.09186
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.0054661, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.00813253
## sample estimates:
## outlier frequency (expected: 0.0015863453815261 )
## 0
testDispersion(so.ls5_sr, alternative = "less") #again, under dispersed
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.0054661, p-value < 2.2e-16
## alternative hypothesis: less
so.ls5a <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = compois)
AIC(so.ls5a) #3290
## [1] 3289.953
so.ls5a_sr <- simulateResiduals(so.ls5a, n = 1000, plot = TRUE)
testResiduals(so.ls5a_sr) #still underdispersed
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.042601, p-value = 0.3267
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.03035, p-value = 0.002
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.01109438
## sample estimates:
## outlier frequency (expected: 0.00152610441767068 )
## 0
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.042601, p-value = 0.3267
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.03035, p-value = 0.002
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.01109438
## sample estimates:
## outlier frequency (expected: 0.00152610441767068 )
## 0
so.ls5b <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = genpois)
AIC(so.ls5b) #3274.5
## [1] 3274.506
so.ls5b_sr <- simulateResiduals(so.ls5b, n = 1000, plot = TRUE)
## qu = 0.75, log(sigma) = -2.866691 : outer Newton did not converge fully.
testResiduals(so.ls5b_sr) #still underdispersed
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.030648, p-value = 0.7378
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.11734, p-value = 0.002
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.01109438
## sample estimates:
## outlier frequency (expected: 0.0015863453815261 )
## 0
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.030648, p-value = 0.7378
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.11734, p-value = 0.002
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.01109438
## sample estimates:
## outlier frequency (expected: 0.0015863453815261 )
## 0
Test second model
so.ls6 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(so.ls6) #3380
## [1] 3380.052
so.ls6_sr <- simulateResiduals(so.ls6, n = 1000, plot = TRUE)
testResiduals(so.ls6_sr) # fails dispersion
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.052567, p-value = 0.1275
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.0083402, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 0.92
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.00813253
## sample estimates:
## outlier frequency (expected: 0.00178714859437751 )
## 0
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.052567, p-value = 0.1275
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.0083402, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 0.92
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.00813253
## sample estimates:
## outlier frequency (expected: 0.00178714859437751 )
## 0
so.ls6a <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = genpois)
AIC(so.ls6a) #3275.5
## [1] 3275.479
#compois distro fails and won't produce DHARMa results
so.ls6a_sr <- simulateResiduals(so.ls6a, n = 1000, plot = T)
testResiduals(so.ls6a_sr) # fails dispersion
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.03162, p-value = 0.7019
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.11438, p-value = 0.002
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.01004016
## sample estimates:
## outlier frequency (expected: 0.00160642570281124 )
## 0
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.03162, p-value = 0.7019
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.11438, p-value = 0.002
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.01004016
## sample estimates:
## outlier frequency (expected: 0.00160642570281124 )
## 0
so.ls6 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(so.ls6)
## [1] 3380.052
so.ls6a <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = genpois)
AIC(so.ls6a)
## [1] 3275.479
# compois distro fails, won't produce DHARMa results
so.ls6b <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = tweedie)
AIC(so.ls6b)
## [1] 3246.35
so.ls6b_sr <- simulateResiduals(so.ls6b, n = 1000, plot = TRUE)
testDispersion(so.ls6b_sr) #fails
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.096452, p-value = 0.006
## alternative hypothesis: two.sided
so.ls6c <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls6c)
## [1] 3877.962
so.ls6c_sr <- simulateResiduals(so.ls6c, n = 1000, plot = TRUE)
testResiduals(so.ls6c_sr)
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.031202, p-value = 0.7174
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1, p-value = 1
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 1, observations = 498, p-value = 0.6306
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 5.083768e-05 1.113682e-02
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.002008032
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.031202, p-value = 0.7174
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1, p-value = 1
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 1, observations = 498, p-value = 0.6306
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 5.083768e-05 1.113682e-02
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.002008032
testDispersion(so.ls6c_sr) #passes
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1, p-value = 1
## alternative hypothesis: two.sided
testQuantiles(so.ls6c_sr)
##
## Test for location of quantiles via qgam
##
## data: simulationOutput
## p-value = 0.4331
## alternative hypothesis: both
so.ls6d <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = gaussian)
AIC(so.ls6d)
## [1] 3943.262
so.ls6d_sr <- simulateResiduals(so.ls6d, n = 1000, plot = TRUE) #fails
plot(so.ls6d_sr)
testResiduals(so.ls6d_sr)
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.095627, p-value = 0.0002216
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0131, p-value = 0.858
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 5, observations = 498, p-value = 0.003539
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.003267824 0.023273861
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.01004016
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.095627, p-value = 0.0002216
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0131, p-value = 0.858
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 5, observations = 498, p-value = 0.003539
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.003267824 0.023273861
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.01004016
so.ls1 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls1) # 3880.6
## [1] 3880.595
# test piri ba
so.ls2 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls2) #3879.6
## [1] 3879.555
lrtest(so.ls1, so.ls2) #p = 0.3
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.BA_piri +
## l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1 |
## Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral +
## avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 16 -1924.3
## 2 15 -1924.8 -1 0.9601 0.3272
# test total BA
so.ls3 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls3) #3885.4
## [1] 3885.409
lrtest(so.ls1, so.ls3) # p = 0.009
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.BA_piri +
## l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1 |
## Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_piri + l.Mineral +
## avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 16 -1924.3
## 2 15 -1927.7 -1 6.8144 0.009043 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# test mineral
so.ls4 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls4) #3881.5
## [1] 3881.533
lrtest(so.ls2, so.ls4) # p = 0.046
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral +
## avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + avgLD_l +
## l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 15 -1924.8
## 2 14 -1926.8 -1 3.9781 0.0461 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# test avg ld
so.ls5 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls5) #3877.7---- keep avg LD, lower AIC
## [1] 3877.655
lrtest(so.ls5, so.ls2) # p = 0.48
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral +
## l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral +
## avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 14 -1924.8
## 2 15 -1924.8 1 0.1004 0.7514
# test piri
so.ls6 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls6) #3882.6
## [1] 3882.59
lrtest(so.ls6, so.ls4) # p = 0.08
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + avgLD_l + l.Veg_Total +
## offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + avgLD_l +
## l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 13 -1928.3
## 2 14 -1926.8 1 3.0575 0.08036 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# test other
so.ls7 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.BA_HA + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls7) #3891
## [1] 3891.009
lrtest(so.ls6, so.ls7) # p = 0.001
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + avgLD_l + l.Veg_Total +
## offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.BA_HA + avgLD_l + l.Veg_Total + offset(l.TFT) +
## (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 13 -1928.3
## 2 12 -1933.5 -1 10.418 0.001248 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# test veg
so.ls8 <- glmmTMB(Shrub_Oak ~ Treat_Type +l.other + l.BA_HA + avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; false convergence (8). See vignette('troubleshooting'),
## help('diagnose')
#won't converge
lrtest(so.ls6, so.ls1) #p = 0.046
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + avgLD_l + l.Veg_Total +
## offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.BA_piri +
## l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1 |
## Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 13 -1928.3
## 2 16 -1924.3 3 7.9957 0.0461 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#issue is that I'm representing 3 measurement scales. reduce BA?
so.ls6 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
ls.all3$BA_10m <- (ls.all3$BA_HA/1000)
ls.all3$l.BA_10m <- log(ls.all3$BA_10m + 1)
so.ls9 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_10m + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; false convergence (8). See vignette('troubleshooting'),
## help('diagnose')
so.ls10 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls10) #3884.5
## [1] 3884.51
so.ls11 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls11) #3395.4
## [1] 3895.37
lrtest(so.ls10, so.ls11) #keep veg
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.other + avgLD_l + l.Veg_Total + offset(l.TFT) +
## (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.other + avgLD_l + offset(l.TFT) +
## (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 12 -1930.2
## 2 11 -1936.7 -1 12.86 0.0003357 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
so.ls12 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + avgLD_l + l.BA_10m + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls12) #3883.8
## [1] 3883.846
lrtest(so.ls11, so.ls12) #keep ba
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.other + avgLD_l + offset(l.TFT) +
## (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.other + avgLD_l + l.BA_10m + offset(l.TFT) +
## (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 11 -1936.7
## 2 12 -1929.9 1 13.523 0.0002357 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#maybe if i drop something else model will run?
so.ls13 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.BA_10m + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls12) #3883.8
## [1] 3883.846
so.ls1a <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_10m + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls1a) # 3873.8
## [1] 3873.765
# this is killing me - moving the basal area down to the 10m scale, to keep measurements on 2 scales; overwriting models
so.ls1 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_10m + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls1) # 3873.8
## [1] 3873.765
so.ls2 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_10m + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls2)
## [1] 3873.773
so.ls3 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_10m + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls3) #3875.7 without mineral
## [1] 3875.745
so.ls3b <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_10m + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls3b) #3871.9 w/o LD
## [1] 3871.895
lrtest(so.ls2, so.ls3) # p = 0.46
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_10m + l.Mineral +
## avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_10m + avgLD_l +
## l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 15 -1921.9
## 2 14 -1923.9 -1 3.9722 0.04626 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lrtest(so.ls2, so.ls3b) #drop ld for sure
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_10m + l.Mineral +
## avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_10m + l.Mineral +
## l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 15 -1921.9
## 2 14 -1922.0 -1 0.1213 0.7277
# test piri
so.ls4 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_10m + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls4) #3871.9
## [1] 3871.997
#test other
so.ls5 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.BA_10m + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls5) #3881.5
## [1] 3881.538
lrtest(so.ls4, so.ls5) # keep other
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.other + l.BA_10m + l.Mineral + l.Veg_Total +
## offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.BA_10m + l.Mineral + l.Veg_Total +
## offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 13 -1923.0
## 2 12 -1928.8 -1 11.541 0.0006809 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# test veg
so.ls6 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_10m + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls6) #3878
## [1] 3877.987
so.ls7 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls7) #3881.7, keep ba
## [1] 3881.709
so.ls4 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_10m + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls4)
## [1] 3871.997
so.ls4_sr <- simulateResiduals(so.ls4, n = 1000, plot = TRUE)
testResiduals(so.ls4_sr)
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.031202, p-value = 0.7174
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1, p-value = 1
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 1, observations = 498, p-value = 0.6306
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 5.083768e-05 1.113682e-02
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.002008032
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.031202, p-value = 0.7174
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1, p-value = 1
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 1, observations = 498, p-value = 0.6306
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 5.083768e-05 1.113682e-02
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.002008032
plot(so.ls4_sr)
hist(ls.all3$Shrub_Oak)
ggplot(ls.all3, aes(x=Shrub_Oak))+
geom_histogram(binwidth = 2)+
facet_grid(rows = vars(Treat_Type))